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The full-density-matrix numerical renormalization group has evolved as a systematic and transpar- 
ent setting for the calculation of thermodynamical quantities at arbitrary temperatures within the 
NRG framework. It directly evaluates the relevant Lehmann representations based on the complete 
basis sets introduced by Anders and Schiller 1 (2005). In addition, specific attention is given to the 
possible feedback from low energy physics to high energies by the explicit and careful construction 
of the full thermal density matrix, naturally generated over a distribution of energy shells. Specific 
examples are given in terms of spectral functions (fdmNRG), time-dependent NRG (tdmNRG), Fermi- 
Golden rule calculations (fgrNRG), as well as the calculation of plain thermodynamic expectation 
values. Furthermore, based on the very fact that, by its iterative nature, the NRG eigenstates are 
naturally described in terms of matrix product states, the language of tensor networks has proven 
enormously convenient in the description of the underlying algorithmic procedures. This paper 
therefore also provides a detailed introduction and discussion of the prototypical NRG calculations 
in terms of their corresponding tensor networks. 
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I. INTRODUCTION 

The numerical renormalization group (NRG) 2 4 is the 
method of choice for quantum impurity models. These 
consist of an interacting local system coupled to non- 
interacting typically fermionic baths, which in their com- 
bination can give rise to strongly correlated quantum- 
many-body effects. Through its renormalization group 
(RG) ansatz, its collective finite size spectra provide a 
concise snapshot of the physics of a given model from 
large to smaller energies on a logarithmic scale. A rich 
set of NRG analysis is based on these finite size spec- 
tra, including statistical quantities that can be efficiently 
computed within a single shell approach at an essentially 
discrete set of temperatures tied to a certain energy shell. 
Dynamical quantities such as spectral functions, however, 
necessarily require to combine data from all energy scales. 
Since all NRG iterations contribute to a single final curve, 
traditionally it had not been clear how to achieve this in 
a systematic clean way, specifically so for finite temper- 
atures. 

The calculation of spectral properties within the NRG 
started with Oliveira and Wilkins 5 ' 6 (1981, 1985) in the 
context of X-ray absorption spectra. This was extended 
to spectral functions at zero temperature by Sakai et al. 7 
(1989). Finite temperature together with transport prop- 
erties, finally, was introduced by Costi and Hewson 8 
(1992). An occasionally crucial feedback from small to 
large energy scales finally was taken care of by the explicit 
incorporation of the reduced density matrix for the re- 
mainder of the Wilson chain (DM-NRG) by Hofstetter 9 
(2000). While these methods necessarily combined data 
from all NRG iterations to cover the full spectral range, 
they did so through heuristic patching schemes. More- 
over, in the case of finite temperature, these methods had 
been formulated in a single-shell setup that associates a 



well-chosen characteristic temperature that corresponds 
to the energy scale of this shell. For a more complete 
listing of references, see Bulla et al. 4 (2008). 

The possible importance of a true multi-shell frame- 
work for out-of-equilibrium situations had already been 
pointed out by Costi 10 (1997). As it turns out, this can 
be implemented in a transparent systematic way using 
the complete basis sets, which where introduced by An- 
ders and Schiller 1 (2005) for the feat of real-time evolu- 
tion within the NRG (TD-NRG). This milestone devel- 
opment allowed for the first time to use the quasi-exact 
method of NRG to perform real-time evolution to expo- 
nentially long time scales. It emerged together with other 
approaches to real-time evolution of quantum many-body 
systems such as the DMRG. 11,12 While more traditional 
single-shell formulations of the NRG still exist for the 
calculation of dynamical quantities using complete ba- 
sis sets, 1,13 the latter, however, turned out significantly 
more versatile. 14 18 In particular, a clean multi-shell for- 
mulation can be obtained using the full-density-matrix 
(FDM) approach to spectral functions fdmNRG. 14 As it 
turns out, this essentially generalizes the DM-NRG 9 to 
a clean black-box algorithm, with the additional bene- 
fit that it allows to treat arbitrary finite temperatures 
on a completely generic footing. Importantly, the FDM 
approach can be easily adapted to related dynamical cal- 
culations, such as the time dependent NRG (tdmNRG), or 
Fermi- Golden-rule calculations (fgrNRG), 19 21 as will be 
described in detail in this paper. 

For the FDM approach, the underlying matrix prod- 
uct state (MPS) structure of the NRG 14 ' 22 provides an 
extremely convenient framework. It allows for an effi- 
cient description of the necessary iterative contractions 
of larger tensor networks, i.e. summation over shared 
index spaces. 23 Moreover, since this quickly can lead in 
complex mathematical expressions if spelled out explic- 



2 



itly in detail, it has proven much more convenient to 
use a graphical representation for the resulting tensor 
networks. 23 In this paper, this is dubbed MPS diagram- 
matics. It compactly describes the relevant procedures 
that need to be performed in the actual numerical simu- 
lation, and as such also represents a central part of this 
paper. 

The paper then is organized as follows: the remainder 
of this section gives a brief introduction to the NRG, com- 
plete basis sets, its implication for the FDM approach, 
and the corresponding MPS description. As a corollary, 
this section also discusses the intrinsic relation of energy 
scale separation, efficiency of MPS, and area laws. Sec- 
tion II gives a brief introduction to MPS diagrammatics, 
and its implications for the NRG. Section III provides 
a detailed description of the FDM algorithms fdmNRG, 
tdmNRG, as well as fgrNRG in terms of their MPS dia- 
grams. This also includes further additional corollaries, 
such as the generic calculation of thermal expectation 
values. Section IV provides summary and outlook. A 
short appendix, finally, comments on the treatment of 
fermionic signs within tensor networks, considering that 
NRG typically deals with fermionic systems. 



A. Numerical renormalization group and quantum 
impurity systems 

The generic quantum impurity system (QIS) is de- 
scribed by the Hamiltonian 

#QIS = H[ mp + H cp \({f 0fI }) + if bath? (1) 

which consists of a small quantum system (the quantum 
impurity) that is coupled to a non-interacting macro- 
scopic reservoir H hath = Efc M ^ M 4 M 2 fcM' e ' g ' a Fermi 
sea. Here cj, creates a particle in the bath at energy e kfl 
with flavor /i, such as spin or channel, and energy index 
k. Typically, e kfl = s k . The state of the bath at the loca- 
tion r = of the impurity is given by f 0fl = ^^Z k V k c kfI 
with proper normalization A/" 2 = ^Z k V k . The coeffi- 
cients Vk are determined by the hybridization coefficients 
of the impurity as specified in the Hamiltonian [e.g. see 
Eq. (4b) below]. The coupling H C p\({fofi}) then can act 
arbitrarily within the impurity system, while it interacts 
with the baths only through f^J , i. e. its degrees of free- 
dom at the location of the impurity. Overall, the Hilbert 
space of the typically interacting local Hamiltonian Ho 
in Eq. (1) is considered small enough so it can be easily 
treated exactly numerically. 

The presence of interaction enforces the treatment 
of the full exponentially large Hilbert space. Within 
the NRG, this consists of a systematic state-space dec- 
imation procedure based on energy scale separation, 
(i) The continuum of states in the bath is coarse 
grained relative to the Fermi energy using the discretiza- 



tion parameter A > 1, such that with W the half- 
bandwidth of a Fermi sea, this defines a set of intervals 
±jy[A-( m - 2+1 )/ 2 ,A-( m -^/ 2 ], each of which is eventu- 
ally described by a single fermionic degree of freedom 
only. Here m is a positive integer, with the additional 
constant z G [0, 1[ introducing an arbitrary shift, 24 ' 25 to 
be referred to as z-shift. (ii) For each individual flavor 
/x then, the coarse grained bath can be mapped exactly 
onto a semi-infinite chain, with the first site described 
by /o/x and exponentially decaying hopping amplitudes 
t n along the chain. This one-dimensional linear setup is 
called the Wilson chain, 2 

N 

H N = ffo + EE(*»-iti,.L + H - c ')' ( 2 ) 

/i n=l 

where Hqis — lini/v-Kx> Hn- For larger n, it quickly 
holds 15 ' 25 

u n = lim t n _x = ^gp^ ^A-f, (3) 

where uj n describes the smallest energy scale of a Wilson 
chain including all sites up to and including site n (de- 
scribed by f nfJL ) for arbitrary A and z-shift. In practice, 
all energies at iteration n are rescaled by the energy scale 
uj n and shifted relative to the ground state energy of that 
iteration. This is referred to as rescaled energies. 

From the point of view of the impurity, the effects of 
the bath are fully captured by the hybridization function 
Y(e) = 7rp(s)V 2 (s), which is assumed spin-independent. 
For simplicity, a flat hybridization function is assumed 
throughout, i.e. T(e) = T$(W — \e\), with the discretiza- 
tion following the prescription of Zitko and Pruschke 25 . 
If not indicated otherwise, all energies are specified in 
units of the (half-) bandwidth, which implies W := 1. 

1. Single impurity Anderson model 

The prototypical quantum impurity model applica- 
ble to the NRG is the single impurity Anderson model 
(SI AM). 26 29 It consists of a single interacting fermionic 
level (d- level), i.e. the impurity, 

#im P = ^2 £d °™ d ° + U ^d^n di . (4a) 

(7 

with level-position Sdo- an d onsite interaction U. This 
impurity is coupled through the hybridization 

ffcp! = £ (4E v ^ + H - c -) ( 4b ) 

cr k 

V v ' 

= V 7T " /0ct 

to a single spinful non-interacting fermi sea, with T 
the total hybridization strength. Here o\ (c^ ka ) cre- 
ates an electron with spin a G {T? 4-} a ^ the d-level (in 
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the bath with energy index fc), respectively. Moreover, 
n da = d\d a , and h ka = c\ a c k(j . At average occupa- 
tion with a single electron, the model has three physical 
parameter regimes that can be accessed by tuning tem- 
perature: the free orbital regime (FO) at large energies 
allows all states at the impurity from empty to doubly 
occupied, the local moment regime (LM) at intermediate 
energies with a single electron at the impurity and the 
empty and double occupied state at high energy only ac- 
cessible through virtual transitions, and the low-energy 
strong coupling (SC) fixed-point or Kondo regime, where 
the local moment is fully screened by the electrons in the 
bath into a quantum-many-body singlet. 

B. Complete basis sets 

Within the NRG, a complete many-body basis 1 can be 
constructed from the state space of the iteratively com- 
puted NRG eigenstates H n \s) n = E%\s) n . With the NRG 
stopped at some final length N of the Wilson chain, the 
NRG eigenstates w.r.t. to site n < N can be comple- 
mented by the complete state space of the rest of the 
chain, |e) n , describing sites n+1, . . . , N. The latter space 
will be referred to as the environment, which due to en- 
ergy scale separation will only weakly affect the states 
\s) n . The combined states, 

\se) n \s) n <g) |e) n , (5) 

then span the full Wilson chain. Within the validity of 
energy scale separation, one obtains 1 

H N \se) n ~E?\se) n , (6a) 

i.e. the NRG eigenstates at iteration n < N are, to a 
good approximation, also eigenstates of the full Wilson 
chain. This holds for a reasonably large discretization 
parameter A > 1.7. 2 ' 4 ' 30 

With focus on the iteratively discarded state space, this 
allows to build a complete many-body eigenbasis of the 
full Hamiltonian, 1 

l(^ N ) = J2 \se)™ (se\, (6b) 

se,n 

where d^d N describes the full many-body Hilbert space 
dimension of the Hamiltonian H^. Here d refers to the 
state space dimension of a single Wilson site, while do 
refers to the state space dimension of the local Hamil- 
tonian Ho, which in addition to fo also fully incorpo- 
rates the impurity [cf. Eq. (1)]. It is further assumed, 
that the local Hamiltonian Ho is never truncated, i.e. 
truncation sets in for some n = no > 0. Therefore, by 
construction, the iterations n' < no do not contribute 
to Eq. (6b). At the last iteration n = TV, all states are 
considered discarded by definition. 1 The truncation at in- 
termediate iterations, finally, can be chosen either w.r.t. 
to some threshold number Nk of states to keep, while 
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Figure 1: Iterative construction of complete basis set 1 within 
the NRG by collecting the discarded state spaces \s)n from all 
iterations n < N (black space at the left of the gray blocks). 
For a given iteration n, these are complimented by the envi- 
ronment \e) n for the rest of the system n > n, i.e. starting 
from site n + 1 up to the overall chain length N considered 
(gray blocks). In a hand- waving picture, by adding site n+1 
to the system of sites n < n, this site introduces a new lowest 
energy scale to the system, with the effect that existing levels 
become split within a narrow energy window (indicated by the 
spread of levels from one iteration to the next). The impurity, 
and also the first few sites can be considered exactly with a 
manageable total dimension of its Hilbert space still. Yet as 
the state space grows exponentially, truncation quickly sets 
in. The discarded state spaces then, when collected, form a 
complete basis. At the last iteration, where NRG is stopped, 
by definition, all states are considered discarded. 



nevertheless respecting degenerate subspace, or, prefer- 
entially, w.r.t. to an energy threshold Ek in rescaled 
energies [cf. Eq. (3)]. The latter is a dynamical scheme 
which allows for a varying number of states depending 
on the underlying physics. 

The completeness of the state space in Eq. (6b) can 
be easily motivated by realizing that at every NRG trun- 
cation step, by construction, the discarded space (eigen- 
states at iteration n with largest energies) is orthogonal 
to the kept space (eigenstates with lowest eigenenergies) . 
The subsequent refinement of the kept space at later iter- 
ations will not change the fact, that the discarded states 
at iteration n remain orthogonal to the state space gener- 
ated at later iterations. This systematic iterative trunca- 
tion of Hilbert space while building up a complimentary 
complete orthogonal state space is a defining property of 
the NRG, and as such depicted schematically in Fig. 1. 



C. Identities 

This section deals with notation and identities related 
to the complete basis sets within the NRG. These are 
essential when directly dealing with Lehmann represen- 



4 



tations for the computation of t her mo dynamical quanti- 
ties. While the combination of two basis sets discussed 
next simply follows Ref. 1, this section also introduces 
the required notation. The subsequent Sec. ID then de- 
rives the straightforward generalization to multiple sums 
over Wilson shells. 

Given the complete basis in Eq. (6), it holds 1 



N 



(7) 



Here the state space projectors are defined to project 
into the kept (X = K) or discarded (X = D) space of 
Wilson shell n. This then allows to rewrite Eq. (7) more 
compactly as 



pK = 



(AT) 

E^ 

n'>n 



(8) 



where the upper limit in the summation, n' < N is im- 
plied if not explicitly indicated. With this, two indepen- 
dent sums over Wilson shells can be reduced into a single 
sum over shells, 1 



E 

m ,U2 



pD f>D _ pDpD , pDpD i V" p D p D 

r n\ ri2 / j n n / v n i n / v n n<i 

(n 1 =n 2 )=n n 1 >(n 2 =n) (ni=n)<n 2 



pKpD 

n ^ n 



5D pK\ 
n 1 n J 



(9) 




For simplified notation, the prime in the last single sum 
over Wilson shells (Yin) indicates that also the kept- 
sectors are included in the sum over Wilson shells, yet 
excluding the all- kept sector XX' ^ KK, since this sector 
is refined still in later iterations. 1,14 

While Eq. (9) holds for the entire Wilson chain, exactly 
the same line of arguments can be repeated starting from 
some arbitrary but fixed reference shell n, leading to 



pKpK 



N 

E 

n\ ,n2>n 



N 



pD pD _ V^' pXpX' 

ri\ r i'2 / j h h ' 



(10) 



where Eq. (8) was used in the first equality. 



D. Generalization to multiple sums over shells 

Consider the evaluation of some physical correlator 
that requires m > 2 insertions of the identity in Eq. (6b) 
in order to obtain a simple Lehmann representation. Ex- 
amples in that respect are tdmNRG or (higher-order) cor- 
relation functions, as discussed later in the paper. In 



all cases, the resulting independent sum over arbitrarily 
many identities as in Eq. (6b) can always be rewritten 
as a single sum over Wilson shells. The latter is desir- 
able since energy differences, such as they occur in the 
Lehmann representation for correlation functions, should 
be computed within the same shell, where both con- 
tributing eigenstates are described with comparable en- 
ergy resolution. 

Claim: Given m full sums as in Eq. (6b), this can be 
rewritten in terms of a single sum over a Wilson shell n, 
such that Eq. (9) generalizes to 

N N ^Ki...K m 

E 4 D 1 i -.-4 D r = E E ''I' ■■■>'»"■ (id 

ni,...,n m h Xf-X m 



where again the prime in the last single sum over Wilson 
shells (Yin) m dicates that all states are to be included 
within a given iteration n, while only excluding the all- 
kept sector Xi, . . . , X m ^ K, . . . , K. Note that via Eq. 
(8), the l.h.s. of Eq. (11) can be rewritten as 

pKi pK m = pD pD 

^no-l " • n -l / j rii ' ' rim 



where no > is the first iteration where truncation sets 
in. This way, Pn -i refers to the full Hilbert space still. 
Proving Eq. (11) hence is again equivalent to proving for 
general n that 



pK m = pD pD = pX 1 AX m 

ni ,...,n m >n n>n 



(12) 

with the upper limit for each sum over shells, rii < N 
and ft < TV, implied, as usual. Therefore the sum in the 
center term, for example, denotes an independent sum 

Yn x >n f ° r aU U i Wltn * = 1j • • • ' m - 

Proof: The case of two sums (m = 2) was already 
shown in Eq. (10). Hence one may proceed via induction. 
Assume, Eq. (12) holds for m — 1. Then for the case m, 
one has in complete analogy to Eq. (9), 

pKi pK m _i . pK m = 
n ' ' n n 

= (E'^^i(E^) 



n'>n 

/ _j n n 

h>n 



UP? 

V n 



pKi pKm-ipDf, 



p^ 1 pXm 

/ j r h • • • r h ' 



- E ^ 

h>n 

where from the second to the third line, it was used that, 

E'E = E' + E' + E' . 

n'>n n m >n n<(n=n'=nm) n<(n=n')<nm n<(fi=n m )<n / 

and the last term in the third line followed from the in- 
ductive hypothesis. This proves Eq. (12). ■ 
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Alternatively, the m independent sums over 
{ni,...n m } in Eq. (12) can be rearranged such, 
that for a specific iteration n, either one of the indices 
rti may carry ft as minimal value, while all other sums 
range from > n. This way, by construction, the index 
rii stays within the discarded state space, while all other 
sums fly are unconstrained up to > = n, thus 
represent either discarded at iteration ft or discarded at 
any later iteration which corresponds to the kept space 
at iteration n. From this, Eq. (12) also immediately 
follows. 



E. Energy scale separation and area laws 

By construction, the iterative procedure of the 
NRG generates an MPS representation for its energy 
eigenbasis. 22 ' 31 This provides a direct link to the den- 
sity matrix renormalization group (DMRG), 32,33 and 
consequently also to its related concepts of quantum 
information. 23 For example, it can be demonstrated that 
quite similar to the DMRG, the NRG truncation w.r.t. 
to a fixed energy threshold E^ is also quasi-variational 
w.r.t. to the ground state of the semi-infinite Wilson 
chain. 15,30 Note furthermore that while DMRG typically 
targets a single global state, namely the ground state of 
the full system, at an intermediate step nevertheless it 
also must deal with large effective state spaces describ- 
ing disconnected parts of the system. This again is very 
much similar to the NRG, which at every iteration needs 
to deal with many states. 

Now the success of variational MPS, i.e. DMRG, to 
ground state calculations of quasi-one-dimensional sys- 
tems is firmly rooted in the so-called area law for the 
entanglement or block entropy Sa = tr(— pA log pa) with 
pA = tr#(/?). 34 36 In particular, the block entropy Sa 
represents the entanglement of some contiguous region 
A with the rest B of the entire system A U B consid- 
ered. This allows to explain, why MPS, indeed, is ideally 
suited to efficiently capture ground state properties for 
quasi-one-dimensional systems. 

In constrast to DMRG for real-space lattices, however, 
NRG references all energy scales through its iterative di- 
agonalization scheme. It zooms in towards the low energy 
scales ( "ground state properties'''' ) of the full semi-infinite 
Wilson chain. Therefore given a Wilson chain of suf- 
ficient length TV, without restricting the case, one may 
consider the fully mixed density matrix built from the 
ground state space |0)jv of the last iteration, for simplic- 
ity. This then allows to analyze the entanglement entropy 
S n of the states \s) n , i.e. the block of sites n' < n, with 
respect to its environment \e) n . The interesting conse- 
quence in terms of area law is that one expects the (close 
to) lowest entanglement entropy S n for the stable low- 
energy fixed point, while one expects S n to increase for 
higher energies, i.e. with decreasing Wilson shell index 
n. 

This is nicely confirmed in a sample calculation for the 
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Figure 2: NRG and area law - analysis for the symmetric 
SI AM for the parameters as shown in panels (b) and (c) [cf. 
Eq. (4); all energies in units of bandwidth]. Panel (a) shows 
the standard energy flow diagram of the NRG for even iter- 
ations where the different colors indicate different symmetry 
sectors. Panel (b) shows the entanglement entropy S n of the 
Wilson chain up to and including site n < N with the rest of 
the chain, given the overall ground state (N = 99). Due to 
intrinsic even-odd alternations, even and odd iterations n are 
plotted separately. Panel (c) shows the actual number of mul- 
tiplets kept from one iteration to the next, using a dynamical 
truncation criteria w.r.t. a predefined fixed energy threshold 
Ek as specified. The calculation used SU(2) sp i n ® SU(2) c har g e 
symmetry, hence the actual number of kept states is by about 
an order of magnitude larger [e.g. as indicated with the max- 
imum number of multiplets kept, Nk in panel (c): the value 
in brackets gives the corresponding number of states] . 



SIAM, as demonstrated in Fig. 2. Figure. 2(a) shows the 
standard NRG energy flow diagram (collected finite size 
spectra, here for even iterations), which clearly outlines 
the physical regimes of free orbital (FO, n < 25), local 
moment (LM, 25 < n < 60), and strong coupling (SC, 
n > 60) regime. Here in order to have a sufficiently wide 
FO regime, a very small onsite interaction U was cho- 
sen relative to the bandwidth of the Fermi sea. Panel 
(b) shows the entanglement entropy S n between system 
(n f < n) and environment (p! > n). Up to the very be- 
ginning or the very end of the actual chain (the latter is 
not shown), this shows a smooth monotonously decaying 
behavior vs. energy scale. In particular, consistent with 
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the area law for lowest-energy states, the entanglement is 
smallest once the stable low-energy fixed point is reached. 
Having chosen a dynamical (quasi- variational) 15 trunca- 
tion scheme w.r.t. to a threshold energy Ek in rescaled 
energies [cf. Eq. (3)], the qualitative behavior of the en- 
tanglement entropy is also reflected in the number of 
states that one has to keep for some fixed overall ac- 
curacy, as shown in Fig. 2(c). Clearly, up to the very 
few first shells prior to truncation, the largest number 
of states must be kept at early iterations. While this 
is a hand- waving argument, this nevertheless confirms 
the empirical fact, that the first few Wilson shells with 
truncation are usually the most important, i.e. most ex- 
pensive ones. Therefore for good overall accuracy all the 
way down to the low energy sector, one must allow for 
a sufficiently large number of states to be kept at early 
iterations. 

The entanglement entropy as introduced above to- 
gether with the area law thus is consistent with the en- 
ergy scale separation along the Wilson chain in [cf. Fig- 
ure. 2(b)]. However, note that the specific value of the 
entanglement entropy is not a physical quantity, in that 
it depends on the discretization. While the entangle- 
ment entropy clearly converges to a specific value when 
including a sufficient number of states, it nevertheless 
sensitively depends on A. The smaller A, the larger the 
entanglement entropy S n is going to be, since after all, 
the Wilson chain represents a gapless system. The over- 
all qualitative behavior, however, is expected to remain 
the same, i.e. independent of A. Similar arguments 
hold for entanglement spectra and their corresponding 
entanglement flow diagram, which provide significantly 
more detailed information still about the reduced den- 
sity matrices constructed by the bipartition into system 
and environment. 15 



F. Full density matrix 

Given the complete NRG energy eigenbasis |se)^, the 
full density matrix (FDM) at arbitrary temperature T = 
1//3 is simply given by 14 

^FDM (r) = ^^? |se) DD H) (13) 

sen 

with Z(T) = J^ne sgd e ~^ E * • By construction of a ther- 
mal density matrix, all energies E™ from all shells n ap- 
pear on an equal footing relative to a single global energy 
reference. Hence any prior iterative rescaling or shifting 
of the energies E™, which is a common procedure within 
the NRG [cf. Eq. (3)], clearly must be undone. From a 
numerical point of view, typically the ground state energy 
at the last iteration n = N for a given NRG run is taken 
as energy reference. In particular, this ensures numeri- 
cal stability in that all Boltzmann weights are smaller or 
equal 1. 

Note that the energies E™ are considered independent 
of the environmental index e. As a consequence, this 



leads to exponentially large degeneracies in energy for the 
states \se) n . The latter must be properly taken care of 
within FDM, as it contains information from all shells. 
By already tracing out the environment for each shell, 
this leads to 14 

~FDM (T) = £^£^|«< S |, (14) 

n s 

y v * v ' 

=w n =P^{T) 

with d the state-space dimension of a single Wilson site, 
and the proper normalization by the site-resolved par- 
tition function Z n (T) = ^2 seDn e~P E * of the density 
matrices p^(T) built from the discarded space of a spe- 
cific shell n only. Therefore, tr(p^(T)) = 1, and also 
Z ( T ) = T,n Z niT). Equation (14) then defines the 
weights w n , which themselves represent a normalized dis- 
tribution, i.e. ^2 n w n = 1. Importantly, Eq. (14) demon- 
strates that FDM is constructed in a well-defined man- 
ner from a distribution of density matrices p^(T) from a 
range of energy shells n. 

1. Weight distribution w n 

The qualitative behavior of the weights w n can be un- 
derstood straightforwardly. With the typical energy scale 
of shell n given by 

uj n = aK~^\ (15) 

with a some constant of order 1. [cf. Eq. (3)], this allows 
to estimate the weights w n as follows, 

\n(w n ) ~ \n(d N - n e- (3uJn /Z) 

— (N — n) hi{d) — j3uo n + const. 

For a given temperature T, the shell n with maximum 
weight is determined by, 

£ lnK) ~ - ln(d) + A -n/2 I Qj 

with the solution 

since the second term is 1//3 times some constant of order 
1. This shows that the weight distribution w n is strongly 
peaked around the energy scale of given temperature T. 
With T = aA -nT / 2 and therefore n T — n*, the distribu- 
tion decays super-exponentially fast towards larger en- 
ergy scales n«nT (dominated by e~^ n with exponen- 
tially increasing uo n with decreasing n). Towards smaller 
energy scales n > n^, on the other hand, the distribu- 
tion w n decays in a plain exponential fashion (dominated 
by d~ n , since with f3u n <C 1, e~^ n 1). 

An actual NRG simulation based on the SIAM is shown 
in Fig. 3. It clearly supports all of the above qualitative 
analysis. It follows for a typical discretization parameter 
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Figure 3: Typical FDM weight distribution calculated for 
the SIAM [cf. Eq. (4)] for the parameters as shown in the 
panel and temperature T = 10 -6 (all energies in units of 
bandwidth). The maximum number of states Nk kept at ev- 
ery iteration was taken constant. The distribution is strongly 
peaked around the energy shell n* > nr, where nr (indicated 
by vertical dashed line) corresponds to the energy scale of 
temperature as defined in the text. The inset plots the weights 
w n on a logarithmic scale, which demonstrates the generic 
plain exponential decay for small energies n > nr, and super- 
exponentially fast decay towards large energies (n < nr)- 



A and local dimension that tit is slightly smaller than 
n*, i.e. towards larger energies to the left of the maxi- 
mum in w ni typically at the left onset of the distribution 
w n , as is seen in the main panel in Fig. 3 (tit is indi- 
cated by the vertical dashed line). Within the shell n* of 
maximum contribution to the FDM, therefore the actual 
temperature is somewhat larger relative to the energy 
scale of that iteration [note that this is clearly related 
to the factor /3, 3 ' 4 introduced by Krishna-murthy et al. 3 
(1980) on heuristic grounds for the optimal discrete tem- 
perature representative for a single energy shell]. 

An important practical consequence of the exponential 
decay of the weights w n for n ^> tlt is that by taking a 
long enough Wilson chain to start with, fdmNRG automat- 
ically truncates the length of the Wilson chain at several 
iterations past ut- Therefore the actual length of the 
Wilson chain N included in a calculation should be such 
that the full distribution w n is sampled, which implies 
that w n has dropped again at least down to ^ 10 -3 . 

The weights w n are fully determined within an NRG 
calculation, and clearly depend on the specific physical 
as well as numerical parameters. Most obviously, this 
includes the state space dimension d of a given Wilson 
site, and the discretization parameter A. However, the 
weights w n also sensitively depend on the specific num- 
ber of states kept from one iteration to the next. For 
example, the weights are clearly zero for iterations where 
no truncation takes place, which is typically the case for 
the very first NRG iterations which include the impurity. 
However, the weights also adjust automatically to the 



specific truncation scheme adopted, such as the quasi- 
variational truncation based on an energy threshold E^. 
In the case of fixed Nk =512 as in Fig. 3, note that if 
d = 4 times the number of states had been kept, i.e. 
Nk = 512 — >• 2048, this essentially would have shifted 
the entire weight distribution in Fig. 3 by one iteration to 
the right to lower energy scales, resulting in an improved 
spectral resolution for frequencies uj < T. 14 For the latter 
purpose, however, it is sufficient to use an increased Nk 
at late iterations only, where around the energy scale of 
temperature the weights w n contribute mostly. 

Furthermore, given a constant number A^k of kept 
states in Fig. 3, the weights w n show a remarkably smooth 
behavior, irrespective of even or odd iteration n. This 
is somewhat surprising at first glance, considering that 
NRG typically does show pronounced even-odd behavior. 
For example, for the SIAM (see also Fig. 2), at even iter- 
ations an overall non-degenerate singlet can be formed to 
represent the ground state. Having no unpaired spin in 
the system, this typically lowers the energy more strongly 
as compared to odd iterations which do have an unpaired 
spin. Therefore while even iterations show a stronger en- 
ergy reduction in its low energy states, its ground state 
space consists of a single state. In contrast, for odd it- 
erations the energy reduction by adding the new site is 
weaker, yet the ground state space is degenerate, assum- 
ing no magnetic field (Kramers degeneracy). In terms of 
the corresponding weight distribution for the full density 
matrix then, both effects balance each other, such that 
distribution of the FDM weights w n results in a smooth 
function of the iteration n, as seen in Fig. 3. 

In summary, above analysis shows that the density ma- 
trix generated by FDM is dominated by several shells 
around the energy scale of temperature. The physical 
information encoded in these shells can critically affect 
physical observables at much larger energies. This con- 
struction therefore shall not be shortcut in terms of the 
density matrix in the kept space at much earlier itera- 
tions, i.e. by using H\s)^ c± Eg\s)^ with the Boltz- 
mann weights thus determined by the energies of the 
kept states. This can fail for exactly the reasons already 
discussed in detail with the DM-NRG construction by 
Hofstetter 9 (2000): the low-energy physics can have im- 
portant feedback to larger energy scales. To be specific, 
the physics at the low-energy scales on the order of tem- 
perature can play a decisive role on the decay channels 
of high-energy excitations. As a result, for example, the 
low-energy physics can lead to a significant redistribution 
of spectral weight in the local density of states at large 
energies. 



2. FDM representation 

The full thermal density matrix p^P M in Eq. (14) rep- 
resents a regular operator with an intrinsic internal sum 
over Wilson shells. When evaluating thermodynamical 
expressions then, as seen through the discussions Sec. I C 
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and ID, its matrix elements must be calculated both with 
respect to discarded as well as kept states. While the for- 
mer are trivial, the latter require some more attention. 
All of this, however, can be written compactly in terms 
of the projections in Eq. (7). 

The reduced density matrix p^P M is a scalar operator, 
from which it follows, 



pX -FDM pX' _ c p^< 



(17) 



This defines the projections R% of p^P M onto the space 
X G {K, D} at iteration n, which are not necessarily nor- 
malized hence the altered notation. Like any scalar op- 
erator, thus also the projections R^ carry a single label 
X only. The projection into the discarded space, 



(18) 



by construction, is a fully diagonal operator as defined 
in Eq. (14). In kept space, however, the originally diag- 
onal FDM acquires non- diagonal matrix elements, thus 
leading to a non-diagonal scalar operator, 



irL = 



pK -FDM p t 



w n ,pKp»(T)pK, 



(19) 



n'>n 



i (T) 



with the properly normalized reduced density matrices, 
[ (T) = tr (/&(r)). (20) 



-FDM/ 
Pn.n' 



These are defined for n' > n and, w.r.t. to the basis 
of iteration n, are fully described within its kept space. 
Note that in the definition of the pP/(T) in Eq. (14) 
the environment consisting of all sites ft > n' had al- 
ready been traced out, hence in Eq. (20) only the sites 
n = n + 1, . . . , n' remain to be considered. By definition, 
the reduced density matrices (T) are built from the 
effective basis |s)°, at iteration n', where subsequently 
the local state spaces of sites n = n' ,n' — 1, . . . , n + 1 
are traced out in an iterative fashion. 

The projected FDM operators R ni like other operators, 
are understood as operators in the basis \s) n , i.e. R* = 
XlsGx(^n)ss / \ s )nn(s'\ (note the hat on the operator), 
while the bare matrix elements (R^) ss f = n (s\R^\s f ) n 
are represented by R% (by convention, written without 
hats). Overall then, the operator R n can be written in 
terms of two contributions, (i) the contribution from iter- 
ation n' = n itself (encoded in discarded space), and (ii) 
the contributions of all later iterations n' > n (encoded 
in kept space at iteration n), 



Rn 



n'>n 



(21a) 



(21b) 



In the last equation, for simplicity, the definition of p n ^ 
for n' > n in Eq. (20) has been extended to include the 
case n' = n, where p n?n = pP(T). 



II. MPS DIAGRAMMATIC S 

Given the complete basis sets which, to a good ap- 
proximation, are also eigenstates of the full Hamilto- 
nian, this allows to evaluate correlation functions in a 
text-book like fashion based on their Lehmann represen- 
tation. Despite the exponential growth of the many- 
body Hilbert space with system size, repeated sums 
over the entire Hilbert space nevertheless can be eval- 
uated efficiently, in practice, due to the one-dimensional 
structure of the underlying MPS [the situation is com- 
pletely analogous to the product, say, of N matrices , 
n G {1, . . . , N}, of dimension D, (A^A™ . . . i^)^ = 

id) There 



Z^fci = l 2-^k 2 = l ' ' ' 2^k N = l ^iM kx,k 2 ' ' ' k N -!,j- 

the sum over intermediate index spaces fci, . . . , fejv-i m 
principle also grows exponentially with the number of 
matrices. By performing the matrix product sequentially, 
however, this is no problem whatsoever]. 



A. Basics and conventions 

The NRG is based on an iterative scheme: given an (ef- 
fective) many-body eigenbasis |s) n _i up to and includ- 
ing site n — 1 on the Wilson chain, a new site with a 
d- dimensional state space \o) n is added. Exact diagonal- 
ization of the combined system leads to the new eigen- 
states 



W)n\s) r 



(22) 



n'>n 



Here the coefficient space As^} lySn of the underlying uni- 
tary transformation is already written in standard MPS 
notation. 23,33 It will be referred to as A-tensor A n which, 
by construction, is of rank 3. Equation (22) is depicted 
graphically in Fig. 4(a): two input spaces (s n _i and a n 
to the left and at the bottom, respectively), and one out- 
put space s n , as indicated by the arrows. Since by con- 
vention in this paper, NRG always proceeds from left to 
right, A-tensors always have the same directed structure. 
Therefore, for simplicity, all arrows will be skipped later 
in the paper. Furthermore, the block A n , which depicts 
the coefficients of the A-tensor at given iteration, will be 
shrunk to a ternary node, resulting in the simplified ele- 
mentary building block for MPS diagrams as depicted in 
panel Fig. 4(b). Finally, note that the start of the Wilson 
chain does not represent any specific specialization. The 
effective state space from the previous iteration is simply 
the vacuum state, as denoted by the (terminating) thick 
dot at the left of Fig. 4(c). The vacuum state represents 
a perfectly well-defined and normalized state, such that 



9 



(a) 

l s )n-l- 



(d) 



An 

T 

\v)n 



(g) 



Sn-1 



(b) 

s) n *n-l . 

s> (e) 
* (h) 

Sn-1 



(c) 
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Figure 4: Basic MPS diagrammatics - Panel (a) Iteration 
step in terms of A-tensor. The coefficient space (A-tensor) for 
given iteration n is denoted by A n , and incoming and outgo- 
ing state spaces are indicated by arrows. Panel (b) cleaned 
up simplified version of diagram in panel (a). Panel (c) in- 
dicates the first A-tensor, in case it has the vacuum state to 
its left, which is denoted by a (terminating) thick dot. Here, 
trivially \o)d = \s)d [with \s)o for n — generated in the 
very next iteration with a Wilson chain in mind]. Panels (d) 
demonstrates the ort ho normality condition of an A-tensor, 
T,aJA [(Tn] y A [an] = 1 [cf. Eq. (23)]. Panel (e) again is fully 
equivalent to panel (d). Panel (f) depicts a reduced density 
matrix. Panel (g) represents the evaluation of matrix ele- 
ments of a local operator B at site n in the effective state 
space s n . Panel (h) again is a cleaned up simplified version of 
panel (g). Panel (i) is similar to panels (g-h), except that the 
operator B was assumed to act at earlier sites on the Wilson 
chain, such that here B already describes the matrix elements 
in the effective basis s n -i, and hence contracts from the left. 
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Figure 5: Basic MPS diagrammatics in the presence of non- 
abelian symmetries - Panel (a) Representation of an irre- 
ducible operator B that acts within the local basis a n in the 
effective basis s n . Being an irreducible operator, a third open 
index emerges, both for the representation of the local op- 
erator B (right incoming index to B an ^ a ' n ^ = (cr n |5 M |cr^)), 
as well as for the overall contracted effective representation 
with the open indizes 5^^^^, where \i identifies the spinor 
component in the irreducible operator B. Panel (b) Sim- 
plified version of panel (a), but exactly the same otherwise. 
Panel (c) Contraction into a scalar representation of an op- 
erator B in the effective representation s n -i which acts at 
some site n < n with operator which acts at site n. With 
B - W = ^ Bp • Bp & scalar operator, the result is a scalar 
operator of rank two in the indices (s n , s' n ). 



all subsequent contractions in the remainder of the pan- 
els in Fig. 4 apply identically without any specific further 
modification. 

Panel Fig. 4(d) depicts the elementary contraction that 
represents the orthonormality condition, 



°s n ,s / n 



i(s\s') n 

£ 4 



«' si*. 



(23) 



again, with panel (e) a cleaned up version, but other- 
wise exactly the same as panel (d). By graphical con- 
vention, contractions, i.e. summation over shared index 
or state spaces, are depicted by lines connecting two ten- 
sors. Note that in order to preserve the directedness of 
lines in Fig. 4(d), it is important w.r.t. bra-states, that 
all arrows on the A* -tensor belonging to bra-states are 
fully reversed. For the remainder of the paper, however, 
this is of no further importance. 

The contraction in panels Fig. 4(d+e) therefore results 
in an identity matrix, given that all input spaces of the 
A-tensor are contracted. For a mixed contraction, such 
as one input and one output state space, on the other 
hand, as indicated in Fig. 4(f), this results in a reduced 
density matrix. There the sum over the state space s n 
is typically weighted by some normalized, e.g. thermal, 
weight distribution p s , as indicated by the short dash 
across the line representing s n together with the corre- 
sponding weights p s . 

Panel Fig. 4(g-i) describe matrix representations of an 
operator B in the combined effective basis s n for a local 
operator acting within a n (panels g-h), or for an opera- 
tor that acted at some earlier site, such that it already 
exists in the matrix representation of the basis s n _i. For 
the latter case, the contraction in panel (h) typically oc- 
curred at some earlier iteration, with subsequent itera- 
tive propagation of the matrix elements as in panel (i) 
for each later iteration. Contracting first the the oper- 
ator B as represented in the state space of a r n [s' n _^\ in 
Fig. 4(h[i]), respectively, followed by the simultaneous 
contraction of (s n _i,cr n ), the cost of the contraction in 
panels (g-i) scales like (D(D 3 ), where D represents the 
matrix dimension for the state spaces s n _i and s n (here 
considered to be the same, for simplicity). 

For the NRG it is crucially important to use 
abelian and non-abelian symmetries for numerical 
efficiency. 2 ' 17,18 ' 37 39 Figure 5 therefore presents elemen- 
tary tensor contractions in the presence of non-abelian 
symmetries. 37 There the basis transformations in terms 
of the A-tensors A n respect the underlying fusion rules 
for non-abelian symmetries. Moreover, elementary oper- 
ators B typically become irreducible operator sets {B^} 
which are described in terms of a spinor with operator 
components labeled by the index /i [see Fig. 5(a)]. Using 
Wigner-Eckart theorem, the arrows, for example, with 
the operator B in Fig. 5(a) imply the underlying Clebsch- 
Gordan coefficient (<j n |/i, <j^). 37 In case of a scalar oper- 
ator 5, the spinor reduces to a single operator, hence \i 
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reduces to a singleton dimension which can be stripped. 
In that case, the third index to the center right of panels 
Fig. 5(a-b) can be removed, resulting in the equivalent di- 
agrams in Fig. 4(g-h). Panel Fig. 5(c), finally, shows the 
contraction of two irreducible operator sets into a scalar 
operator B • £t = £ • B\ , again represented in the 

combined effective basis s n . Here the operator B is con- 
sidered to have acted once at some earlier site, whereas 
its daggered version acts on the current local site n. Note 
that again the daggered (conjugated) version has all its 
arrows reversed where, in addition, in the MPS diagram 
the dagger indicates, that the operator B^ as compared 
to B has already been also flipped upside down. 

The only essential difference when using non-abelian 
symmetries with MPS diagrammatics is the emergence 
of extra indices (lines) w.r.t. to irreducible operators (in- 
dex \i above). The underlying A-tensors, of course, need 
to respect the fusion rules of the symmetries employed, 
but on the level of an MPS diagram, this is implied. A 
detailed introduction to non-abelian symmetries and its 
application to the NRG has been presented in Ref. 37. 
Therefore for the rest of this paper, for simplicity, no fur- 
ther reference to non-abelian symmetries will be made, 
with all tensor networks based on the elementary con- 
tractions already presented in Fig. 4. 

III. FDM APPLICATIONS 

A. Spectral functions 

Consider the retarded Green's function 

G* ( t ) = -i#(t)(B(t)&) T (24) 

=G BC (t) 

which may regarded as the first term in the standard 
fermionic Green's function G R (t) = -i$(t)({B(t)C^}) T . 
Here B(t) = e lHt Be~ lHt , where, as usual, the Hamilto- 
nian H of the system is considered time- independent. In 
Eq. (24), an operator & acts at time t = on a sys- 
tem in thermal equilibrium at temperature T, described 
by the thermal density matrix p(T) = e~^ H /Z(T) with 
(•) T = tr(p(T) • ). The system then evolves to some time 

t > 0, where a possibly different operator B is applied. 
The overlap with the original time evolved wave func- 
tion then defines the retarded correlation function of the 
two events. Fourier-transformed into frequency space, 
G R (u) = / ^e iuJt G R (t), its spectral function is defined 
by 

A bc (uj) ee -IlmGg c M = / £e**G% c (t) 

= J §e iu;t tT^p(T)e iAt Be- iAt C^. (25) 

When evaluated in the full many-body eigenbasis, in 
principle, this requires the insertion of two identities, (i) 



to evaluate the trace, and (ii) in between the operators 
B and & to deal with the exponentiated Hamiltonian. 
For simplified, with the eigenbasis sets 1 = ^ a l a )( a l = 
| b) (b |, the spectral function becomes, 

A B c(oo) = E/ £e i{u - Eah)t Pa{a\B\b)(b\&\a) 

ab 

= Y,PaBabC* ab -S(uj-E ab ), (26) 

ab 

with ee E b — E a and p a = ^e~^ Ea . By convention, as 
usual, operators carry hats, while matrix representations 
in a given basis have no hats (B vs. B a b). Equation (26) 
is referred to as the Lehmann representation of the corre- 
lation function in Eq. (24). In the case of equal operators, 
B = C, the spectral function is a strictly positive func- 
tion, i.e. a spectral density. In either case, the integrated 
spectral function results in the plain thermodynamic ex- 
pectation values, 

J dojA BC (oj) = Y,Pa B ab C* ab = <BC t ) T . (27) 

ab 

Now using the complete NRG eigenbasis, \a) — ^ \se) n 
and | b) — >• |sV) n /, one may have been tempted of di- 
rectly reducing the double sum in Eq. (26) to a sin- 
gle sum over Wilson shells using Eq. (9). This im- 
plies that the thermal weight would be constructed as 
paiT) - e-P Ea e _/3£; "' X from both, the discarded 
(X = K) as well as the kept (X = K) space at itera- 
tion n. This, however, ignores a possible feedback from 
small to large energy scales which has been shown to be 
crucial in the NRG context. 9 

The solution is to take the FDM as it stands in 
Eq. (13). This, however, introduces yet another inde- 
pendent sum c over Wilson shells, in addition to a and b 
in Eq. (26) above, 

A BC (u) = Y,PcaB ab C: c -5(uj-E ab ). (28) 

abc 

The triple-sum over {a, 6, c} can be treated as in Eq. (11). 
With {a,6,c} {s,s',s p } n e {XX% ^ KKK}, nev- 
ertheless, X = X p are locked to each other since p it- 
self represents a scalar operator, i.e. does not mix kept 
with discarded states. Therefore only the contributions 
XX 7 ^ KK as known from a double sum remain. With 

tr($ DM B(t) • Ct) 

= £ E tr(4 x -^ DM .4 x %)Pf 

n XX'X„ V ' 

it follows for spectral functions (fdmNRG), 

AbcH = [Cl Rn] g , s (B n ) ss , 5 (lu- El,), (29) 

n,ss' 
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Figure 6: (Color online) MPS diagram for calculating spec- 
tral functions using fdmNRG based on the Lehmann represen- 
tation in Eq. (29). In general, spectral functions are preceded 
by an NRG forward sweep, which generates the NRG eigen- 
basis decomposition (horizontal lines; cf. Fig. 4). Correlation 
functions then require the evaluation of the matrix elements 
tr(Cp T • ( s )f?( s /)), as indicated at the left of the figure. The 
energies of the indices (states) s and s' are "probed" such that 
their difference determines the energy uj — E™ s , = E™, — E™ 
of an individual contribution to the spectral function, as indi- 
cated by the xS(u — E ss /) next to the indices s and s' in the 
upper right of the figure. The sum X} n '>n m the discarded 
state space of p FDM (T), indicated to the lower right, results 
in the object R n [cf. Eq. (21b)]. The individual contributions 
p™^f (T) are generated by the Boltzmann weights in the dis- 
carded space at iteration n , as indicated to the right. The 
contribution at n — n, i.e. R^, can simply be determined 
when needed. On the other hand, the cumulative contribu- 
tions n > n are obtained in a simple prior backward sweep, 
starting from the last Wilson shell N included, as indicated 
by the small arrow pointing to the left. Having n > n, this 
calculation always maps to the kept space, thus resulting in 
Rn • Finally, the spectral data is collected in a single forward 
sweep, as indicated at the bottom of the figure. 



where the prime with the sum again indicates that only 
the combinations of states ss' G XX' ^ KK are to be 
considered at iteration n. To be specific, given the scalar 
nature of the projections R n , the first term implies the 
matrix product (ClR n ) x ' x = (Cl) x ' x R x . 

The MPS diagram of the underlying tensor structure 
is shown in Fig. 6. Every leg of the "ladders" in Fig. 6 
corresponds to an NRG eigenstate (MPS) \s) n for some 
intermediate iteration n. The blocks for the MPS coef- 
ficient spaces (A-tensors) are no longer drawn, for sim- 
plicity [cf. Fig. 4]. The outer sum over the states s f 
in Eq. (29) corresponds to the overall trace. Hence the 
upper- and lower-most leg in Fig. 6 at iteration n carry 
the same state label s f , as they are connected by a line 
(contraction). Furthermore, the inserted identity in the 
index s initially also would have been identified with two 
legs [similar to what is seen in Fig. 7 later]. At iteration 
n, however, the state space s directly hits the FDM, lead- 
ing to the overlap matrix x (s\s) x = 5 S sb~xx [ nence this 
eliminates the second block from the top in Fig. 7] . As a 



result, only the single index s from the second complete 
sum remains in Fig. 6. The same argument applies for 
the index s" . 

The two legs in the center of Fig. 6, finally, stem from 
the insertion of the FDM which can extend to all itera- 
tions n' > n. Note that the case n' < n does not appear, 
since there the discarded state space used for the con- 
struction of the FDM is orthogonal to the state space s 
at iteration n. The trace over the environment at itera- 
tion n leads to the reduced (partial) density matrices R n . 
Here the environmental states \e) n t for the density ma- 
trices Pn>(T) for n' > n had already all been traced out, 
as pointed out with Eq. (14). The FDM thus reduces at 
iteration n to the scalar operator R^ as introduced in 
Eq. (21). 

In summary, by insisting on using the FDM in Eq. (13) 
this only leads to the minor complication that R% needs 
to be constructed and included in the calculation. The 
construction of R^, on the other hand, can be done in 
a simple prior backward sweep, which allows to generate 
R% iteratively and thus efficiently. All of the R% need 
to be stored for the later calculation of the correlation 
function. Living in kept space, however, the computa- 
tional overhead is negligible. The actual spectral data, 
finally, is collected in a single forward sweep, as indicated 
in Fig. 6. 



1. Exactly conserved sum rules 

By construction, FDM allows to exactly obey sum- 
rules for spectral functions as a direct consequence 
of Eq. (27) and fundamental quantum mechanical 
commutator relations. For example, after complet- 
ing the Green's function in Eq. (24) to a proper 
many-body correlation function for fermions, Gd(t) = 
— (F})t, with cv creating an electron in level 
d at the impurity and {•, •} the anticommutator, the in- 
tegrated spectral function results in 



dujA(uj) = ({dJ^}) T = 1, 



(30) 



due to the fundamental fermionic anticommutator rela- 
tion, {d,(F} = 1. In practice, Eq. (30) is obeyed ex- 
actly within numerical double precision noise (10 -16 ), 
which underlines the fact that the full exponentially large 
quantum-many body state space can be dealt with in 
practice, indeed. Note, however, that Eq. (30) holds by 
construction, and therefore it is not measure for conver- 
gence of an NRG calculation. The latter must be checked 
independently. 15 



2. Implications for complex Hamiltonians 

The Hamiltonians analyzed by NRG are usually time- 
reversal invariant, and therefore can be computed using 
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non-complex numbers. In case the Hamiltonian is not 
time-reversal invariant, i.e. the calculation becomes in- 
trinsically complex, the A-tensors on the lower leg of the 
ladders for the operators B and in Fig. 6 must be 
complex conjugated (see also Fig. 4). Consequently, this 
implies for the FDM projections i? n , that in Fig. 6 its 
constituting A-tensors in the upper leg need to be com- 
plex conjugated. 

B. Thermal expectation values 

Arbitrary thermodynamic expectation values can be 
calculated within the fdmNRG framework, in priniciple, 
through Eq. (27). Given the spectral data on the l.h.s. 
of Eq. (27), for example, this can be integrated to obtain 
the thermodynamic expectation value on the r.h.s. of 
Eq. (27). In practice, this corresponds to a simple sum 
of the non-broadened discrete spectral data as obtained 
from fdmNRG. Using the plain discrete data has the ad- 
vantage that it does not depend on any further details of 
broadening procedures which typically would introduce 
somewhat larger error bars otherwise. 

For dynamical properties within the NRG, usually only 
local operators are of interest. That is, for example, the 
operators B or C in Eq. (27) act within the local Hamil- 
tonian Ho [Wilson shell n = 0; cf. Eq. (1)], or within 
the very first Wilson sites n < no, where no stands for 
the first Wilson shell where truncation sets in. For this 
early part of the Wilson chain, the weights w n are iden- 
tically zero. Consequently, the reduced thermal density 
matrix is fully described for iteration n < no for arbi- 
trary temperatures T by R^(T) in kept space. For a 
given temperature, the aforementioned simple backward 
sweep to calculate R% then already provides all necessary 
information for the simple evaluation of the thermal ex- 
pectation value of any local operator C [e.g. C := BC^ 
inEq. (27)], 

{C) T = tr[i?( K )(T)C( KK >], (n<n ) (31) 

with Cn KK ^ the matrix elements of the operator C in 
the kept space of iteration n. With no truncation yet 
at iteration n, the kept space is the only space avail- 
able, i.e. represents the full state space up to iteration 
n [hence the brackets around the K's]. For strictly local 
operators acting within the state space of i^o, one has 
(C) T = tr[i?o(r)Co]. The clear advantage of Eq. (31) is, 

that once R^(T) has been obtained for given temper- 
ature, any local expectation value can be computed in 
a simple manner without the need to explicitly calculate 
the matrix elements of the operator C throughout the 
entire Wilson chain. 

In Eq. (31), it was assumed that the operator C acts 
on sites n < no only. This can be relaxed significantly, 
however, assuming that temperature is typically much 
smaller than the bandwidth of the system. In that case, 
the weight distribution w n already also has absolutely 



negligible contribution at earlier iterations n' <C tlt 
which clearly stretches beyond no (see Fig. 3 and dis- 
cussion). Hence Eq. (31) can be relaxed to all iterations 
n for which E n '<n w n' < 1- 

In the case that the operator C is not a local opera- 
tor at all, but nevertheless acts locally on some specific 
Wilson site n, then using Eqs. (14) and (21) it follows for 
the general case, 

(C) T = tr[RK(T)C™]+tr[R»(T)C™] 

+ c^%, (32) 

n' <n 

which corresponds to the partitioning of Eq. (14) given by 
= En'>n + En'=n + En'<n> respectively. The last 
term in Eq. (32) derives from the discarded state spaces 
for Wilson shells n' < n at (much) larger energy scales. 
Therefore the fully mixed thermal average applies, such 
that the resulting constant c = ^tr crri ((7) is the plain 

average of the operator C in the local basis \<j n ) that it 
acts upon. To be specific, this derives from the trace over 
the environmental states \e) n in Eq. (14). Equation (31) 
finally follows from Eq. (32), in that for n < no, by con- 
struction, due to the absence of truncation the second 
and third term in Eq. (32) are identically zero. 

C. Time-dependent NRG 

Starting from the thermal equilibrium of some initial 
(I) Hamiltonian H 1 , at time t = a quench at the lo- 
cation of the quantum impurity occurs, with the effect 
that for t > the time-evolution is governed by a differ- 
ent final (F) Hamiltonian H F . While initially introduced 
within the single-shell framework for finite temperature, 1 
the same analysis can also be straightforwardly general- 
ized to the multi-shell approach of fdmNRG. Thus the 
description here will focus on the FDM approach. 

Given a quantum quench, the typical time-dependent 
expectation value of interest is 

Cit) = <C-(t)) T = tr[^(T) • e iHFt Ce- iHFt ], (33) 

with C some observable. While the physically relevant 
time domain concerns the dynamics after the quench, i.e. 
t > 0, one is nevertheless free to extend the definition of 
Eq. (33) also to negative times. The advantage of doing 
so is, that the Fourier transform into frequency space of 
the C(t) in Eq. (33) defined for arbitrary times becomes 
purely real, as will be shown shortly. With this the ac- 
tual time-dependent calculation can be performed in fre- 
quency space first in a simple and for the NRG natural 
way, 

C(uj) = J §e iujt tr(p l (T)-e iAFt Ce- iAFt ) , (34) 

A Fourier transform back into the time-domain at the 
end of the calculation, finally, provides the desired time- 
dependent expectation value C(t) = J C(uj)e~ lUJt duo for 



13 




iteration 



collecting spectral data in a single sweep having (s,s',Sj) g {KKK} 



Figure 7: (Color online) MPS diagram for the calculation of 
quantum quenches using tdmNRG [cf. Eq. (36)]. The calcula- 
tion is performed in frequency space as depicted, which at the 
very end of the calculation is Fourier transformed back into 
the time domain to obtain the desired time-dependent expec- 
tation value C{t) [cf. Eq. (34)]. The calculation requires a 
complete eigenbasis for initial (black horizontal lines) and fi- 
nal Hamiltonian [orange (gray) horizontal lines] , respectively, 
which are computed in two preceding NRG runs. Their re- 
spective shell-dependent overlap matrices S n (two light gray 
boxes at lower left) are calculated in parallel to the calcula- 
tion of the matrix elements C (dark gray box at the top). The 
projections B} n of the FDM (box at the lower right) are evalu- 
ated with respect to the initial Hamiltonian, but have exactly 
the same structure otherwise as already discussed with Fig. 6 
[see also Eq. (21)]. The spectral data, finally, is collected in a 
full forward sweep, as indicated by the arrow at the bottom. 
To be specific, the summation is over all Wilson shells n and 
for a given iteration n, over all states (s, s f , Si) £ {KKK} with 

Si e {si,s 2 }. 

t > 0. In order to obtain smooth data closer to the 
thermodynamic limit, a weak log- Gaussian broadening 
in frequency space quickly eliminates artificial oscilla- 
tions in the time-domain which derive from the loga- 
rithmic discretization. Note that for the sole purpose 
of damping these artificial oscillations, typically a signif- 
icantly smaller log-Gauss broadening parameter a < 0.1 
suffices as compared to what is typically used to obtain 
fully smoothened correlation functions in the frequency 
domain, e.g. a > 0.5 for A = 2 {e.g. see EPAPS of 
Ref. 14). 

1. Lehmann representation 

For the Lehmann representation of Eq. (34), in princi- 
ple, three complete basis sets are required: one completed 
basis set i derived from an NRG run in H l to construct 
p T (T), and two complete basis sets / and /' from an NRG 
run in H F to be inserted right before and after the C op- 
erator, respectively, to describe the dynamical behavior. 
Clearly, two NRG runs in H l and H F are required to de- 



scribe the quantum quench. 5 ' 6 ' 40 With this, the spectral 
data in Eq. (34) becomes 

CM = ^(m&(T)<i\f)C ff ,6(w-E} f ,)(35) 

which generates the overlap matrix S. Now using the 
complete NRG eigenbasis sets together with the FDM, 
again similar to the fdmNRG in Eq. (29), this introduces 
another sum over Wilson shells. Therefore the fourier- 
transformed time-dependent NRG (tdmNRG) becomes, 

n,ss f 

where (s,s') G {X,X'}. In addition, X e {K, D} de- 
scribes the sector of the reduced density matrix R l n from 
the initial system. To be specific, the notation for the 
first term in Eq. (36) implies the matrix product, 

For example, the left daggered overlap matrix S n selects 
the overlap of the sectors {X, X 7 } between initial and final 
eigenbasis, respectively. The prime in Eq. (36) indicates 
that the sum includes all combinations of sectors XX'X ^ 
KKK, i.e. a total of seven contributions. The latter 
derives from the reduction of the independent three-fold 
sum over Wilson shells [cf. Eq. (35)] into a single sum 
over Wilson shells n as discussed with Eq. (11). It is 
emphasized here, that the reduction of multiple sums in 
Wilson shells as in Eqs. (9) and (11) is not constrained 
to having the complete basis sets being identical to each 
other. It is easy to see that it equally applies to the 
current context of different basis sets from initial and 
final Hamiltonian. 

The MPS diagram corresponding to Eq. (36) is shown 
in Fig. 7. It is similar to Fig. 6, yet with several essential 
differences: the block describing the matrix elements of 
the original operator B has now become the block con- 
taining C. The original operator & is absent, i.e. has 
become the identity. Yet since its "matrix elements" are 
calculated with respect to two different basis sets (ini- 
tial and final Hamiltonian), an overlap matrix remains 
(lowest block in Fig. 7). In the context of the correlation 
functions in Fig. 6, the bra-ket states for the inserted 
complete basis set in the index s could be reduced to the 
single bra-index s, such that it affected a single horizontal 
line only. Here, however, two different complete basis sets 
hit upon each other, which inserts another overlap matrix 
(second block from the top in Fig. 7, which corresponds 
to the Hermitian conjugate of the lowest block). The re- 
duced density matrices ^' x , finally, are built from the 
initial Hamiltonian, yet are completely identical in struc- 
ture otherwise to the ones already introduced in Eq. (21). 

The basis of the initial Hamiltonian enters through the 
two legs (horizontal black lines) in Fig. 7 which connect 
to the density matrix R n . All other legs refer to the 
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NRG basis generated by the final Hamiltonian [horizon- 
tal orange (gray) lines]. Finally, note that the plain con- 
traction S^RS of the lower three tensors w.r.t. the in- 
dices si and 52 can simply be evaluated through efficient 
matrix multiplication as in Eq. (37), while nevertheless 
respecting the selection rules on the state space sectors 
{X,X',X}^{K,K,K}. 
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D. Fermi-Golden-Rule calculations 

The NRG is designed for quantum impurity models. 
As such, it is also perfectly suited to deal with local 
quantum events such as absorption or emission of a gen- 
eralized local impurity in contact with non-interacting 
reservoirs. 5,6,16,19,20,40 If the rate of absorption is weak, 
such that the system has sufficient time to equilibrate 
on average, then the resulting absorption spectra are de- 
scribed by Fermi's- Golden rule (fgr), 41 



collecting spectral data in single sweep having (s,s') € {KK} 



Figure 8: (Color online) MPS diagram for the calculation of 
absorption spectra using Fermi's- Golden-rule (fgrNRG) medi- 
ated by the operator & [cf. Eq. (40)]. The two center legs 
(horizontal black lines) refer to the state space of the initial 
Hamiltonian, while the outer legs [horizontal orange (gray) 
lines] refer to the state space of the final Hamiltonian. There- 
fore the matrix elements of & are mixed matrix elements 
between eigenstates of initial and final Hamiltonian. 



A(u>) = 2nJ2pl(T)-\(f\&\i)\ 2 -5(u-E if ), (38) 

where i and / describe complete basis sets for initial and 
final system, respectively. The system starts in the ther- 
mal equilibrium of the initial system. The operator 
describes the absorption event at the impurity system, 
i.e. corresponds to the term in the Hamiltonian that 
couples to the light field. The transition amplitudes be- 
tween initial and final Hamiltonian are fully described 
by the matrix elements Cif = (i\C\f). Given that the 
energy difference Eif = E^ — E\ in Eq. (38) needs to 
be calculated between states of initial and final system, 
absorption or emission spectra usually show threshold 
behavior in the frequency uj. The threshold frequency 
is given by the difference in the ground state energies of 
initial and final Hamiltonian, cj t hr = = Eg — E l g , 

which eventually is blurred by temperature. 

The difference between absorption and emission spec- 
tra is the reversed role of initial and final system, while 
also having & —> C. That is, from the perspective of 
the absorption process, the emission process starts in the 
thermal equilibrium of the final Hamiltonian, with sub- 
sequent transition matrix elements to the initial system. 
This also implies that emission spectra have their contri- 
butions at negative frequencies, i.e. frequencies smaller 
than the threshold frequency uj t hr indicating the emis- 
sion of a photon. Other than that, the calculation of an 
emission spectrum is completely analogous to the calcu- 
lation of an absorption spectrum. With this in mind, 
the following discussion will be therefore constrained to 
absorption spectra only. 

While an absorption spectrum is already defined in 
frequency domain, it nevertheless can be translated into 



the time domain through Fourier transform, 

A{t) = J ^e~^A(u) 

= ^^(T).e^(i|C|/)e-^.(/|(7tK) 

hf 

= (e^ lt Ce-^ Ft -&) l T . (39) 

V v ' 

=C(t) 

Thus absorption spectra can also be interpreted similar 
to correlation functions and quantum quenches: at time 
t = an absorption event occurs (application of the op- 
erator C\ which for example rises an electron from a 
low lying level into some higher level that participates in 
the dynamics). This alters the Hamiltonian, such that 
the subsequent time evolution is governed by the final 
Hamiltonian. At some time t > then, the absorption 
event relaxes back to the original configuration (appli- 
cation of C). Therefore A(t) essentially describes the 
overlap amplitude of the resulting state with the original 
state with no absorption within the thermal equilibrium 
of the initial system. While the "mixed" time evolution 
of C(t) in Eq. (39) may appear somewhat artificial at 
first glance, it can be easily rewritten in terms of a single 
time-independent Hamiltonian: by explicitly including a 
further static degree (e.g. a low lying hole from which 
the electron was lifted through the absorption event, or 
the photon itself), this switches H l to H F , i.e. between 
two dynamically disconnected sectors in Hilbert space of 
the same Hamiltonian (compare discussion of type-1 and 
type-2 quenches in Ref. 21). 

Within the FDM formalism, the calculation of Fermi- 
Golden rule calculations as defined in Eq. (38) becomes 
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(fgrNRG), 

A(oj) = ^'H^n\s' s (Cn)ss'S(cO-E^), (40) 
n,ss f 

where (s, s f ) G {X J ,X F } £ {KK}. Therefore (C n ) M / = 
n^l^kOn represents mixed matrix elements between 
states from initial and final Hamiltonian, respectively, 
which nevertheless can also be easily calculated using the 
basic contractions discussed with Fig. 4. 

The MPS diagram to be evaluated for Eq. (40) is shown 
in Fig. 8. Its structure is completely analogous to the cal- 
culation of generic correlation functions in Fig. 6, except 
that similar to the quantum quench earlier, here again 
the basis sets from two different Hamiltonians come into 
play. 5 ' 6 ' 40 In contrast to the quantum quench situation 
in Fig. 7, however, no explicit overlap matrices are re- 
quired. Instead, all matrix elements of the local oper- 
ator & themselves are mixed matrix elements between 
initial and final system. The reduced density matrices 
R l n are constructed w.r.t. the initial Hamiltonian, but 
again exactly correspond to the ones already introduced 
in Eq. (21) otherwise. 

1. Technical remarks 

Absorption or emission spectra in the presence of An- 
derson orthogonality or strongly correlated low-energy 
physics typically exhibit sharply peaked features close to 
the threshold frequency with clear physical interpreta- 
tion. While in principle, a single Hamiltonian with dy- 
namically disconnected Hilbert space sectors may have 
been used, this is ill-suited for an NRG simulation. Us- 
ing a single NRG run, this can only resolve the low- 
energy of the full Hamiltonian, i.e. of the initial system 
as it is assumed to lie lower in energy. Consequently, 
the sharp features at the threshold frequency will have 
to be smoothened by an energy window comparable to 
co>thr = AE g in order to suppress discretization artifacts. 
This problem is fully circumvented only by using two sep- 
arate NRG runs, one for the initial and one for the final 
Hamiltonian. 5 ' 6 ' 40 With the NRG spectra typically col- 
lected in logarithmically spaced bins, having two NRG 
runs then, it is important that the data is collected in 
terms of the frequencies v = uo — o; t hr taken relative to 
the threshold frequency cj t hr as defined earlier. 



E. Higher-order correlation functions 

Consider the spectral function of a three-point corre- 
lation function, which in the time domain is given by 

Abcd^M) = (D^Cih^T (41) 
ee tr (p(T) • e^De^-^Ce-^B) 
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Figure 9: MPS diagram for the evaluation of a three-point 
correlation functions as in Eq. (42), which thus generalizes 
fdmNRG (see also Fig. 6). 



Given a time-invariant Hamiltonian, the correlator of 
three operators C, and D acting at three different 
times results in the dependence on effectively two times 
t\ and t 2l since to as in B(to) can simply be chosen arbi- 
trary, i.e. to — for simplicity. 

Using the NRG eigenbasis sets together with the FDM, 
the Lehman representation of Eq. (41) requires four in- 
dependent sums over complete eigenbasis sets, one from 
the FDM (X p ), and three by inserting an identity with 
every exponentiated Hamiltonian (Xi, X2, X3 from left to 
right), respectively. Again, with the reduced density ma- 
trix p being a scalar operator, one has X p = Xi. Using 
Eq. (11) then, in frequency space Eq. (41) becomes 

A B Cd(0Ji,0J 2 ) = Yl [ B n R n] s ^ Sl (D n ) Sl s 2 (C n )s 2 s 3 

n S1S2S3 

x S^-E^JS^-E^J, (42) 



with (si, ^2,53) G {Xi,X 2 ,X 3 } 7^ {KKK}, as indicated 
by the prime next to the sum, having [£? n .R n ] 3 ' 1 = 
B^'^Rn 1 • The MPS diagram corresponding to the 
spectral representation in Eq. (42) is shown in Fig. 9. 

The more challenging part with Eq. (42) is the de- 
pendence on two frequencies. While the corresponding 
full collection of data into bins (ui,u 2 ) can become ex- 
pensive, however, certain fixed frequency points together 
with different kernels corresponding to a different ana- 
lytic structure of the higher-order correlation function 
[which then replace the S- functions in Eq. (42)], appear 
feasible with reasonable effort. Moreover, within the 
NRG context, by construction, one has comparable en- 
ergy resolution for uj\ and uj 2 at a given energy shell. 
Hence it remains to be seen in what respect vastly differ- 
ent energy scales of oj\ as compared to uj 2 ^ if required, are 
affected by the ansatz of energy scale separation within 
the NRG. 
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IV. SUMMARY AND CONCLUSIONS 

The framework of tensor network has been applied 
to the NRG. This makes full use of the complete basis 
sets as introduced by Anders and Schiller 1 (2005) which, 
within the approximation of energy scale separation, also 
represent many-body eigenstates of the full Hamilto- 
nian. Together with the full density matrix (FDM) ap- 
proach, complete basis sets allow for simple transpar- 
ent algorithms, as demonstrated for correlation function 
(fdmNRG), time-dependent quenches (tdmNRG), as well as 
Fermi- Golden-rule (fgrNRG) calculations. The underly- 
ing principle is based on the plain Lehmann representa- 
tion of the relevant dynamical expressions, which within 
the NRG, can be evaluated in a text-book like clean and 
transparent fashion. 

The framework of complete basis sets clearly allows 
for straightforward further generalizations. For example, 
one can envisage multiple consecutive time-steps which 
thus generalizes tdmNRG with the possibility to implement 
periodic switching. 42 While initially, the system starts in 
thermal equilibrium of a given Hamiltonian, after each 
quench the description of the system must be projected 
onto the complete basis set of the following Hamiltonian 
in terms of the reduced density matrix. For all these cal- 
culations, however, one must keep in mind that Wilson 
chains are not thermal reservoirs. 43 Within the tdmNRG, 
for example, this can manifest itself as finite size effect, 
in that already for a single quench in the absence of an 
external magnetic field, an initial excess spin at the im- 
purity cannot be fully dissipated into the bath even in the 
limit of time t — »• oo, leading to (small) residual magne- 
tization at the impurity. In cases where these discretiza- 
tion effects become a strong limiting factor, hybrid NRG 
approaches have been devised with the idea to extend the 
bath to a more refined or uniform spectrum. However, 
since this typically compromises energy scale separation 
along the full Wilson chain, other methods such as the 
DMRG need to be incorporated. 22,44 
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Appendix A: Fermionic signs 

The NRG is typically applied to fermionic systems (for 
extensions to bosonic systems see, for example, Refs. [4, 
45,46]). Through its iterative prescription, the resulting 



MPS has a specific natural fermionic order in Fock space, 
\s) n = J2 (A^A^-...-A^) s \a n )...\ao)\a d ), 



= \(Tn,...,(T0,O-d) 



(Al) 

where |cr^) stands for the local state space of the impu- 
rity. Site n' > n is added after site n, hence the state 
space |cr n /) naturally appears to the left \a n ) with sec- 
ond quantization in mind. The environmental states \e) n 
w.r.t. to iteration n which refers to the sites n' > n is 
irrelevant for the following discussion, and hence will be 
skipped. 

Let be a fermionic operator that acts on the impu- 
rity. Here $ is assumed an arbitrary operator that nev- 
ertheless creates or destroys an odd number of fermionic 
particles such that fermionic signs apply. A very fre- 
quent task then is to represent this operator in the effec- 
tive many-body-basis at iteration n, i.e. to calculate the 
matrix elements (C^) ss > = n (s\c^\s f ) n [cf. Fig. 4]. This 
involves the basic matrix-element w.r.t. to local state 
spaces, 

(cr n , ... ,0-0,^1^10-^, ... ,cr f ,cr f d ) = 

= [ n (w-i) n *0] • foi a v*>' ( A2 ) 



,0 s - 



with z = (— l) n = exp(z7rn), akin to the Pauli z-matrix. 
That is, by pulling the operator c\ acting on the impu- 
rity, to the right past the second quantization operators 
that create the states <j n . , fermionic signs apply, resulting 
in a Jordan- Wigner string 



Zi, 



(A3) 



=0,...,' 



to be called z- string in short. Note that through 
the Jordan- Wigner transformation, which maps fermions 
onto spins and vice versa, exactly the same string oper- 
ator as in Eq. (A3) emerges. For a one-dimensional sys- 
tem with nearest neighbor hopping, the Jordan- Wigner 
transformation to spins allows to eliminate on the oper- 
ator level of the Hamiltonian further complications with 
fermionic signs. This is fully equivalent, of course, to 
the explicit treatment of the Jordan- Wigner string in a 
numerical setting that keeps a fermionic basis. The op- 
erators zi in Eq. (A2) take care of the book keeping of 
fermionic signs, by inserting —1 (+1) for all states 
at site i with odd (even) number of particles n (Ji . The 
operators Zi are diagonal and hence commute with each 
other. In the case of additional explicit spin-degrees of 
freedom, such as the localized spin in the Kondo model, 
its z-operator is proportional to the identity matrix and 
hence can be safely ignored. 

In the following, three alternative viewpoints are dis- 
cussed for dealing with fermionic signs in the MPS setup 
of the NRG. To be specific, the following discussion as- 
sumes c* = d) which creates a particle at the impurity's 
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Figure 10: (Color online) MPS diagrams and fermionic signs. 
Consider the matrix elements of a local operator d) which 
creates a particle at the impurity, i.e. the first local state 
space of the MPS in the effective MPS space \s) n . A z-string 
(Jordan- Wigner string) Z = (g) . z\ arises [green (gray) hori- 
zontal line in the middle]. The endpoints (open circles) indi- 
cate the range of the z-string, i.e. starting from and includ- 
ing site to site n. For every crossing of the z-string with a 
black line, which represent state spaces, fermionic signs ap- 
ply. Panel (a) shows that a z-string can be rerouted (light 
dashed lines, pushed in the direction of the red arrow). The 
resulting configuration in panel (b) shows that by rerouting 
the z-string significantly fewer crossings with black lines can 
be achieved. In particular, the z-strings which applied to all 
sites to the right of & (panel a), can be significantly reduced 
to local fermionic signs at the impurity and another fermionic 
sign with the state space s n . 

d-level. As such, it generates a Jordan Wigner string 
for all sites added subsequently to the MPS, i.e. sites 
z = 0,...,n [cf. Eq. (A3)]. 

Viewpoint 1: rerouting of z-string in tensor network 

Figure 10 depicts an MPS diagram for the typical eval- 
uation of matrix elements with relevant fermionic signs. 
The A-tensors that derive from a preceeding iterative 
state space generation of the NRG are depicted by the 
ternary nodes (cf. Fig. 4). By keeping track of the to- 
tal number n of particles for all indices then, for some 
specific index a the fermionic sign is given by ( — l) Ua . 

The z-string that is required for the evaluation of the 
matrix elements of d) , stretches across all local state 
spaces Gi with < i < n. This is depicted by the light 
green (gray) line in Fig. 10 (note that this is not the ex- 
tra index that takes care of non-abelian symmetries as in 
Fig. 5, even though graphically its role is not that dis- 
similar). Here the interpretation is such, that a crossing 
of the z-string with a state space inserts fermionic signs 
for this state space. 47 49 Consider then, for example, the 
upper right A-tensor, A n: in Fig. 10. For simplicity, its 
three legs are labeled / = s n _i (state space from previous 
iteration), a^ n \ (new local state space), and r = s n (com- 



bined state space) for left, local, and right, respectively. 
By tracking the total particle number for all states, given 
the left-to-right orthonormalization (see arrows in Fig. 4), 
by construction it must hold n\ + n a = n r . The index a 
is crossed by the z-string, hence fermionic signs apply at 
the location of the crossing, 

Z G EE (-1)"* = (_l)»r(_l)-n, = ZlZr . (A4) 
= (-l)+»i 

Therefore, instead of applying fermionic signs with in- 
dex a, it is equally correct to apply fermionic signs 
with the indices I and r. This allows to reroute the z- 
string 47 49 as indicated in Fig. 10 (dashed line to the 
upper right, with the shift in the z-string indicated by 
short red arrow). Note that for this rerouting to work, 
the actual left-to-right orthonormalization is not strictly 
required, and could be relaxed, in general, to the more 
general condition n\ ± n r ± n a = even. In particular, 
this includes n\ ± n r ± n a = 0, which suggests that 
any direction of orthonormalization is acceptable, to- 
gether with a generic current site that combines all (ef- 
fective) state spaces to an even number of particles, i.e. 
n\ + n r + na = n to t = even (for n to t = odd, a global 
minus sign would apply in the case of rerouting). 

The basic rerouting step as indicated above can be 
repeated, such that the z-string can be pulled to the top 
outside the MPS diagram in Fig. 10(a), with the final 
configuration shown in Fig. 10(b). The state to the very 
left (black dot) is the vacuum states with no particles, 
hence the z-string can be fully pulled outside at the left. 
As a result, instead of the original n crossings with the 
state space <j n , only two crossings of the z-string with 
state spaces (black lines) remain: one crossing with the 
local state space at the impurity itself, leading to 

$ ->$z d = (zi)i, (A5) 

which fully acts within the state space of the impurity, 
and another crossing with the state space \s f ) n at itera- 
tion n. 

In typical applications which include thermal expecta- 
tion values or correlation functions, however, an operator 
d) never comes by itself, as its expectation value with 
respect to any state with well-defined particle number 
would be zero. Therefore creation and annihilation oper- 
ators always appear in pairs. For the local spectral func- 
tion, for example, d) is paired with its daggered version 
d. In their overall combination, the fermionic signs w.r.t. 
the index s' appear twice and hence annihilate each other. 
This situation is sketched in Fig. 11. The matrix element 
discussed previously with Fig. 10 is shown in the upper 
half of the figure. Given the case of spectral functions 
(cf. Fig. 6), its counterpart is shown at the bottom. The 
reduced density matrix R is a scalar operator, such that 
the particle number of the states s and s" must match. 
Similarly, the outer two states are connected through 
the overall trace (solid line to the very right), hence re- 
fer to exactly the same state. Consequently, the same 
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tors in the full many-body Hilbert space without making 
reference to MPS notation. Given the fermionic order of 
sites as in Eq. (Al), a fermionic operator c& that destroys 
a particle at site k, has the tensor-product form 



lr 



...1 



k-l 



1 c k ® z k +i 



z n , (A6) 



where 1^ is the identity matrix at site i, c k the desired 
operator acting within the state space of site fc, and Zi = 
( — l) n% as in Eq. (A2). Now, applying a z-operator to 
the states s' at the last site n is equivalent to applying a 
z-operator to each individual site, 



Figure 11: (Color online) Example: fermionic signs in cor- 
relation functions. Two MPS diagrams as in Fig. 11 for the 
matrix elements of d and d) are combined, as required, for ex- 
ample, for the calculation of correlation functions. The result- 
ing product of matrix elements n{s'\d\s") ri • R^j s • n {s\(P\s f ) n 
leads to cancelation of the fermionic signs in the index s' in 
the rerouted z-strings (light green lines), as indicated by the 
two splashes to the right. Hence the right end-point of the 
z-string can be fully retracted to the very left of the diagram, 
as indicated by the dashed red arrows. The partial contribu- 
tion R to the FDM is a scalar operator, such that assuming 
charge conservation, the particle number of the states s and 
s" also must be same. Hence the z-string in Fig. 11 could have 
been equally well also rerouted downwards, instead. The re- 
spective fermionic signs with states s and s" still would have 
canceled, while the order of application of the z-operator with 
the impurity would have changed. 



ZCh 



n 

(®*) 



i Ck 



Zd ® zo ® • • • z k -i ® [zc] k <g> i fe +i . . . i n ,(A7) 



since (zi) 2 = 1^. In the application to thermodynamic 
quantities such as correlations functions, the operator C k 
would again appear together with its daggered version 
C\, hence insertion of Z 2 has no effect, yet can be split 
in equal parts, i.e. C\Ck = (ZCkY(ZCk)- Therefore, 
ZCk can be equally well used instead of Ck- As a result, 
similar to Fig. 11, the z-strings have again been fully 
flipped from the sites to the right of site k to the left 
of site k, with the additional transformation Ck —> [zc\k- 
Note that, essentially, this equivalent to fully reverting 
the fermionic order. 



fermionic sign factor applies twice with the rerouted z- 
strings, which thus cancels, i.e. [(— l) ns ] 2 = 1 (indicated 
by the two splashes with s' at the right). Consequently, 
the right end-point of the z-strings can be retracted along 
the rerouted z-string all to the way to the left of the im- 
purity (indicated by the red dashed arrow). 

Therefore given the A-tensors for the basis transfor- 
mations from a prior NRG run that only generates the 
basis, above line of argument allows to ignore fermionic 
signs for most of the subsequent calculation of thermo- 
dynamic quantities or spectral properties. Specifically, in 
given example which applies to fdmNRG, tdmNRG, as well 
as fgrNRG, it is sufficient to calculate the spectral func- 
tions for the operator d — » z&d [cf. Eq. (A5)] and fully 
ignore fermionic signs for the rest of the chain. This is 
in contrast to the original setup where the full z-string 
needs to be included. 



Viewpoint 2: Operator representation 

An alternative way to demonstrate the effect of rerout- 
ing of the z-string can be given by looking at the equiva- 
lent (numerical) tensor-product representation of opera- 



Viewpoint 3: Auxiliary fermionic level 

In the case of absorption spectra, the absorption of a 
photon creates an electron- hole pair, hftoft, where the hole 
hft can be simply treated as a spectator in the dynamics. 
Nevertheless, by explicitly including the hole in the cor- 
relation function, i.e. by using the operator oft — >• Weft, 
this operator itself now already forms a pair of fermions 
that preserves particle number (assuming that hft creates 
a hole). Therefore, by construction, h) oft simply com- 
mutes with all Wilson sites except for the impurity upon 
which it acts. 

The same argument can be repeated for a standard 
spectral function, by introducing an auxiliary fermionic 
level h that does not participate in the dynamics, i.e. 
does not appear in the Hamiltonian. In general, prepend- 
ing the states in Eq. (Al) by the states \ah) of the "hole" , 
i.e. 

\a n , . . . ,a ,a d ) -> \a n , . . . , a , a d )\a h ), (A8) 

immediately results in the same consistent picture as al- 
ready encountered with Fig. 11 or Eq. (A7). 
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